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Abstract. - A model of rod- like polyelectrolyte brushes in the presence of monovalent and multi- 
valent counterions but with no added-salt is studied using Monte Carlo simulation. The average 
height of the brush, the histogram of rod conformations, and the counterion density profile are ob- 
tained for different values of the grafting density of the charge-neutral wall. For a domain of grafting 
densities, the brush height is found to be relatively insensitive to the density due to a competition 
between counterion condensation and inter-rod repulsion. In this regime, multivalent counterions 
collapse the brush in the form of linked clusters. Nematic order emerges at high grafting densities, 
resulting is an abrupt increase of the brush height. 



Electrostatic correlation effects of highly charged macro-ions in aqueous solutions have 
recently attracted much attention [1]. Novel phenomena such as charge inversion and like- 
charge attraction, which cannot be understood within the Poisson-Boltzmann theory, appear 
to play a key role in biological processes such as DNA packaging [2] and cell scaffolding dy- 
namics [3]. Such correlations, and their role in determining the overall structural properties 
of the system, are most pronounced in high density polyelectrolyte solutions where the 
charge "patterns" that form due to the correlations strongly interact with each other. If the 
polyelectrolyte is sufficiently stiff in its backbone structure [such as DNA and filamentous 
actin (F-actin)], a geometric constraint is imposed on the charges in each polyelectrolyte 
segment to adopt a more or less linear configuration, which makes a solution of such rod-like 
polyelectrolytes a patterned matrix for the counterions (that is, the ions of opposite charge 
that are present to neutralize the overall solution). The interplay between these correla- 
tions and patterns can lead to novel structural and dynamical properties [4-7] , the study of 
which could be useful in understanding the biological phenomena mentioned above as well 
as designing functional biomimetic materials. 

A common realization of such high density assemblies occurs in polyelectrolyte brushes 
[8], where charged polymers are end-grafted to surfaces. While polyelectrolyte brushes have 
been mostly studied because of their role in stabilization of colloidal suspensions, under- 
standing their structural properties could help us in a variety of problems in biophysics and 
-technology. A case in point is filamentous microtubules, i.e. rod-like structures build from 
tubulin monomers with charged C-terminal amino-acid tails extending from the filament 
core [9]. The C-terminal tails seem to play key roles in microtubule function, e.g. in con- 
nection with binding of microtubule-associated proteins and processivity of the molecular 
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Fig. 1: Snapshots of the system with monovalent counterions: a) a = 1.2L, b) a = 0.55L and c) 
a — 0.25L. Inset: the schematic configuration of a polyelectrolyte rod of length L and two rotational 
degrees of freedom 9 and <f>. 

motors that move along microtubules [10]. Presumably they form disordered brushes mak- 
ing it difficult to study them by systematic scattering and osmotic stress experiments [11]. 
Design of more effective DNA chips [12] may also benefit from a better understanding of 
rod-like polyelectrolyte brushes. Interestingly STM studies suggest that densely grafted ar- 
rays of short single-stranded DNA molecules on gold substrates self-organize into disordered 
structures, and it takes an applied electric field normal to the grafting surface to align these 
brushes [13]. 

Most theoretical studies of polyelectrolyte brushes have focused on flexible chains, i.e. 
charged polymers that are longer than their persistence lengths. If the polyelectrolytes are 
sufficiently short or highly charged, they would behave effectively as charged rods [14] that 
are constrained to be very close to each other because of the brush structure; a feature 
that can bring about frustration [7]. Here, we study the structural properties of rodlikc 
polyelectrolyte brushes in the presence of monovalent and trivalent counterions using Monte 
Carlo (MC) simulation techniques. We calculate the average thickness of the brush, the 
density profile of the counterions, and the statistics of the conformation of the rods, as 
functions of the surface coverage density. We find three distinct regimes for the collective 
behavior of the rods: (i) a weakly interacting regime at low densities where the counterions 
are dispersed in the solution (see Fig. [T^) , (ii) an intermediate regime where the counterions 
are condensed in the brush (see Fig. Hb), and (iii) a novel nematic regime at high densities 
(see Fig. QJ). The counterion density profile is found to be depleted near the wall in the 
intermediate regime, while the statistics of the rod conformations show that a small fraction 
of the rods choose to lie on the surface as a result of electrostatic frustration. We also find 
that trivalent counterions collapse the brush in the intermediate regime by forming links 
between clusters of neighboring rods (see Fig. [2^). Interestingly, the height of the brush 
in the intermediate regime seems to depends very weakly on the areal density of the rods, 
both for monovalent and trivalent counterions. 

We model each polyelectrolyte as a rigid rod of length L with N m charged hard sphere 
monomers of diameter a and charge q rn = — e distributed along it in a uniform periodic 
array. N r polyelectrolytes are end grafted on a square lattice with a lattice constant a on 
a surface, so that the areal grafting density is p a — I /a 2 . Each polyelectrolyte has two 
rotational degrees of freedom 9 and <f> (Fig. [T^,). The counterions are modeled as hard-core 
spheres of diameter a and they carry a charge q c = Ze in which Z is the valence of the 
counterions. The system is electrostatically neutralized: N c x Z — N r x N m , where N c is 
the number of counterions. 

We use Lekner summation method [15] to account for long-range Coulomb interactions 
by applying lateral periodic boundary conditions in x and y directions (see Fig. QJi), and the 
simulation box is finite in z direction, normal to the brush surface. Elementary moves for 
counterions, are chosen at random in the simulation box around it. For the polyelectrolyte 
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rods the non-grafted end of each rod moves randomly on the surface of a hemisphere of radius 
L. For highly charged rods it is known that algorithms involving independent moves for the 
polyelectrolytcs become inefficient since many counterions are linked to polyelectrolytes, 
and moving a rod away from these counterions would cost a lot of energy that effectively 
prohibits any Monte Carlo move [16]. To solve this problem, we follow Ref. [16] and look 
for all counterions at distances less than d of each rod and rotate them all together with 
the polyelectrolyte. Considering the acceptance rate of the MC moves, the value of d has 
been optimized to d = 2a. We equilibrate the system by starting from different initial 
configurations and reaching the same stationary value of the Coulomb energy. In simulations 
with monovalent counterions, the equilibrium was reached after 2 x 10 MC steps per particle 
for the smallest value of the lattice constant. For the case of trivalent counterions the 
equilibration time is 3 x 10 4 MC steps per particle. The simulation box height in z direction, 
H, has been chosen in such a way that the results are independent of its value; H = 12L 
for monovalent and H = 8L for trivalent counterions. 

To quantify the strength of the electrostatic energy relative to the thermal energy, we 
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can define the Bjerrum length Ib = £fc e BT , where e is the dielectric constant of the solvent. 
In our simulations, we chose to fix the Bjerrum length — which in room temperature is about 
0.7 nm for water — to a fraction of the length of the rod, namely, Zb = L/7. This means 
that the length of rods is chosen to be about L ~ 5 nm at room temperature. We also 
fixed the diameter of the rods to a = 0.1L. In the case of monovalent counterions, the 
simulation parameters were: N m = 7, iV r = 7 x 7 = 49, and N c = 343. Thermal averages 
were obtained over 6 x 10 4 MC steps, after equilibrating the system. In simulations with 
trivalent counterions (Z = 3), we chose: N m = 7, N r = 9 X 9 = 81, and iV c = 189, and 
thermal averages are calculated over 9 x 10 4 MC steps after equilibration. Note that the 
Manning parameter £ = zn ™ 1b [17] was set to £i = 1 for monovalent, and £3 = 3 for 
trivalent counterions. 

We vary the lattice constant a, keeping the other parameters fixed, and calculate the 
brush height as h — jj- J2i=i ( cos #i) > where (• • •} denotes ensemble average. For each value 
of a, we start the system from a random initial configuration and calculate thermal averages 
after equilibration of the system. Figure 3 shows the brush height as a function of the lattice 
constant a for Z = 1. When the lattice constant is very large, the effective attraction of 
the rods to the counterions becomes very weak and the counterions evaporate away from 
them, while the long distance between the rods makes their own electrostatic interaction 
negligible. In this regime, cos 9 has a uniform distribution with the thermal average of i, 
and the brush height is h = L/2. As the spacing a is reduced, the rods start to interact 
with each other electrostatically and the brush is swollen, as can be seen in Fig. [3] Note 
that in this regime the majority of the counterions are still away from the surface (see the 
snapshot of the system in Fig. QJi). To justify this observation, one can roughly consider the 
z = plane as a charged surface with the unform charge density a s — -^JS for which it is 
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Fig. 3: The brush height as a function of the lattice constant a, in the presence of monovalent 
counterions. Inset: the log-log plot of h/L — 0.5 as a function of a for Z = 1. At large lattice 
constants the brush height scales as a~ 2 (the dashed line has the slope of —2). 

Fig. 4: The brush height versus a in the presence of trivalent counterions. 

known from Poisson-Boltzmann theory that the counterions are confined within a distance 

^ = ■2ixZi- B a = iitZIbN ' ca U e d the Gouy-Chapman length. We can readily see that 
A 3> L, so long as a > I and )»a that is a requirement for the smeared surface charge 
density approximation to be justified. 

The swelling of the brush in this regime seems to obey a scaling law of the form \~\ ^ 
^j-, as demonstrated in the inset of Fig. [3] This scaling behavior can be understood by 
assuming that individual rods are fluctuating as dipolcs in the presence of a uniform electric 
field as obtained from the smeared charge density cr s , namely, the sum of monopolc-dipolc 
interactions which are the dominant terms in the multipole expansion valid for large a. 
The scaling behavior is reminiscent of that of flexible polyelectrolyte brushes in the Pincus 
regime [8]. 

Upon decreasing the lattice constant, the effective surface charge density of the brush is 
increased so that the counterions are attracted to the brush, as the snapshot of the system in 
Fig. [TJd shows. This trend continues until a new regime sets in where all the counterions are 
confined within the brush, which is analogous to the osmotic regime in flexible polyelectrolyte 
brushes [8] . One could roughly estimate the onset of this crossover by setting the thickness of 
the counterion layer (~ A) equal to the height of the brush (~ L), which yields a x ~ ^ 1 / 2 L. 
In this regime, the counterions are involved in heavily screening the interaction between the 
rods, and thus one expects that the swelling slows down, as can be seen in the plateau-like 
feature in Fig. [3] for the intermediate range of 0.551/ < a < 1.2L. At higher densities, 
the rods are completely neutralized by the counterions, and the remaining excluded volume 
interaction among them promotes a novel nematic ordering that sets in sharply at a = 0.551/ 
(see Fig. |3|), similar to the Onsager nematic transition in 3D solutions of hard rods [18]. A 
typical snapshot of the system in the nematic phase is shown in Fig. [TJ;. 

The picture becomes more complicated with the introduction of multivalent counterions. 
Figure [4] shows the brush height as a function of the lattice constant in the presence of 
trivalent counterions. For infinitely large a the the brush height is still h = L/2, and it does 
swell initially as a is decreased similar to the monovalent case. However, once a sufficient 
fraction of the trivalent counterions are attracted to the rods and fluctuate in the vicinity 
of them, they give rise to a fluctuation-induced attraction [19]. This attractive interaction 
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Fig. 5: The histogram of cos 8 for a) Z = 1 and a > 1.2L, b) Z = 1 and a < 1.2L, and c) Z = 3. 

competes with the repulsion between the rods and the result of this competition turns 
back the swelling of the brush around a = 2.2L (see Fig. @|. This trend becomes more 
drastic as of a — 2.0L, where the brush height is ft. = Q.5L again, because the trivalent 
counterions can now link the polyelectrolyte rods to each other (see Fig. [2^) and collapse 
the brush [20]. Similar to the monovalent case, an interesting plateau regime appears for 
0.551/ < a < 1.0L, where the enhancement in the repulsion between the rods as a decreases 
is exactly balanced by the increased possibility of forming links between neighboring rods. 
This leads to the clusters growing in size more and more until at a = 0.55L the collection 
of now fully neutralized hard rods undergo a ncmatic transition; the brush height suddenly 
increases (see Fig. [4]) and the rods become predominantly parallel, as the snapshot in Fig. 
[5Jd shows. Note that the onset of the nematic transition coincides in the two cases of mono- 
and trivalent counterions, in line with the suggestion that the transition is entirely driven 
by the excluded volume interaction. 

In addition to the average thickness of the brush, it is instructive to monitor the his- 
togram of the rod conformations P(cos8) in the different regimes. In Fig. [5^i, P(cos6) is 
shown for four different lattice constants in the a > 1.2L domain, for monovalent case. We 
find that in this regime by decreasing a, the cos 8 = part of the histogram decreases and 
the cos 8 = 1 part increases, which corresponds to a "normal" swelling of the brush. (At 
sufficiently large lattice constants the histogram shows a uniform distribution of cos 8, as 
expected.) As the lattice spacing is further reduced to a < 1.2L, the situation becomes 
more subtle, as Fig. [SJd shows. For 0.55L < a < 1.2L the trend of change in the cos (9 = 
part of the histogram is different from Fig. [5^: decreasing the lattice constant a increases 
both the cos 8 = and cos 8 = 1 parts of the histogram, while the intermediate region of 
the histogram decreases to maintain the normalization of the distribution. In this regime, 
the system presumably tries to lower the electrostatic repulsion between the polyelectrolytes 
by increasing angles between the neighboring rods and consequently some of the rods lie 
on the surface, as can be seen in the snapshot for a = 0.55L in Fig. [TJd. For a < 0.55L, 
a drastic change of behavior is observed and the cos 8 = part of the histogram vanishes, 
which signals the nematic ordering of the rods. The histogram is also shown for the case 
of trivalent counterions in Fig. While the histogram is similar to the monovalent case 
in the large-a regime and in the nematic phase, for intermediate values of a it develops a 
maximum at some angle that is characteristic of the clusters that form when the multivalent 
counterions link neighboring rods (see Fig. [2^i) . 

The histogram of cos# can be used to determine the average density profile of the rods. 
One can define the fluctuating density pR d(r) = ^ J2i Jo ds S 3 (r — tjs) where tj de- 
notes the unit director of the i-th rod, and the rods are assumed to have a smeared charge 
distribution for simplicity of the presentation. We define the average density profile as 
n Rod( z ) = jt^i J dxdy (pR d( r )) 5 where the averaging over the direction of the rods can be 
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Fig. 6: The number density of counterions in units of a _2 L _1 as a function of z/L for different 
values of a for a) Z = 1 and b) Z — 3. c) Rescaled number density of monovalent counterions as 
functions of z/X. 

performed using the distribution P(cos#). After some straightforward manipulation, one 
can write the density profile as 7iRod(2) = j^ji f z , L * This expression shows that the 

density profile of the rods will have a singular behavior near the surface if the cos 6 — part 
of the distribution is nonzero, i.e. a non- vanishing fraction of the rods lie on surface. 

We have also obtained the counterion density profile as a function of z for each value 
of a, as shown in Fig. [6^ for Z = 1, and in Fig. [6Jd for Z = 3. For all values of the 
spacing, the general expectation that the counterions are confined within a distance of the 
order the Gouy-Chapman length holds true. This is demonstrated in Fig. [5J; where the 
profiles are compared with the exact solution of the Poisson-Boltzmann equation for a flat 
uniformly charged surface n(z) — 2 -nZ^i B \' 1 (T+zTa) 7 ' with the Gouy-Chapman length being 

2 

A = 2-kZI-bN as a bove. While for a^> L the profiles follow the Poisson-Boltzmann solution, 
deviations appear for a ~ L presumably due to the condensation of the counterions inside 
the brush. It is interesting to note that all the profiles in Fig. collapse on the Poisson- 
Boltzmann solution for z > 5A. Although not attempted here, we note that more elaborate 
solutions of the Poisson-Boltzmann equation suited to the brush geometry can be used to 
account for the profiles obtained via simulation [20]. 

The counterion density profiles appear to develop a depletion in the vicinity of the surface 
for a > 0.55L, as can be seen in Figs. and [HJa. This feature is absent for a < 0.55L, 
presumably due to the emergence of the nematic ordering. In fact, one observes that the 
depletion is present so long as the cos 9 — part of the histogram for the rods is non- 
vanishing. This suggests that the singular behavior of the density profile of the rods near 
the surface (described above) might cause the depletion of the counterions by way of excluded 
volume interaction: as long as some of the rods are lying on the surface, they jam the space 
near the surface and deplete the counterions. The range of the depletion in the profiles of 
Figs. [6^ and[6j3 (~ 2a) is consistent with this picture [21]. 

In conclusion, we have studied the collective behavior of rod-like polyelectrolytes that 
are end-grafted on a substrate, in the presence of mono- and trivalent counterions. We have 
demonstrated that the interplay between condensation of counterions and the repulsion of 
the rods can bring about a number of different regimes, while multivalent counterions can 
introduce correlations that lead to a structured collapse of the brush. The emergence of 
such patterns in the rod conformations — of staggered type for monovalent counterions and 
with "tepee frameworks" in case of the trivalent — seems to originate from similar free energy 
minimizing effects, and might have practical applications, in the design of DNA chips. 

if if if 
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